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We  developed  an  equivalent  capillary  model  of  a  microscale  fiber-fence  structure  to  study  the  microscale 
evolution  and  transport  of  liquid  in  a  porous  media  and  to  reveal  the  basic  principles  of  water  transport  in 
gas  diffusion  layer  (GDL).  Analytical  solutions  using  the  model  show  that  a  positive  hydraulic  pressure  is 
needed  to  drive  the  liquid  water  to  penetrate  through  the  porous  GDL  even  consisting  of  the  hydrophilic 
fibers.  Several  possible  contributions  for  the  water  configuration,  such  as  capillary  pressure,  gravity,  vapor 
condensation,  wettability  and  microstructures  of  the  GDL,  are  discussed  using  the  lattice  Boltzmann 
method  (LBM).  It  is  found  that  the  distribution  manners  of  the  fibers  and  the  spatial  mixed-wettability 
in  the  GDL  also  play  an  important  role  in  the  transport  of  liquid  water. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Owing  to  the  high  efficiency,  low  noise,  and  minimal  pollu¬ 
tion,  the  proton  exchange  membrane  (PEM)  fuel  cell  has  emerged 
as  one  of  the  most  promising  power  sources  for  the  automobile 
engine.  Gas  diffusion  layer  (GDL)  is  one  of  the  essential  compo¬ 
nents  of  the  PEM  fuel  cells.  The  main  roles  of  the  GDL  are  to  allow 
the  gaseous  reactants  (oxygen  in  cathode  and  hydrogen  in  anode) 
to  move  towards  the  catalyst  layer  (CL)  or  the  micro-porous  layer 
(MPL)  and  to  provide  paths  for  liquid  water  to  flow  towards  the 
flow  channel.  It  also  plays  an  important  role  in  the  electron  and 
thermal  transports.  If  a  large  number  of  pores  were  blocked  by 
the  liquid  water  in  the  GDL  the  number  of  the  available  pores  for 
gaseous  flow  would  decrease.  Therefore  the  gaseous  flow  resis¬ 
tance  will  increase  and  the  effective  reaction  area  will  reduce  in 
the  CL.  Control  of  the  liquid  water  distribution  (the  water  manage¬ 
ment  technology),  a  key  technique  to  ensure  high  performance  and 
long  durability  of  PEM  fuel  cells,  has  received  much  attention.  An 
effective  water  management  requires  a  careful  design  of  the  GDL 
microstructure,  such  as  the  porosity,  the  pore  size  distribution,  and 
the  spatial  mixed-wettability  distribution. 

Understanding  the  fundamental  mechanisms  of  water  transport 
in  the  GDL  is  a  key  step  for  the  effective  management  of  liquid 
water.  The  capillary  pressure  and  vapor  condensation  are  consid¬ 
ered  as  two  of  the  main  contributors  to  the  water  configuration. 
Generally,  liquid  water  transport  in  a  porous  media  is  governed  by 
the  capillary  fingering.  It  is  necessary  to  apply  a  pressure  to  force 
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the  liquid  water  to  penetrate  into  the  pores  of  GDL  Therefore,  the 
capillary  pressure  vs.  liquid  water  saturation  curve,  providing  a  con¬ 
stitutive  relationship  to  simulate  the  liquid  water  transport  in  the 
cathode,  has  received  considerable  attention  [1,2].  However,  recent 
studies  show  that  the  vapor  condensation  obviously  affects  the  liq¬ 
uid  water  transport  in  the  cathode  [3-5],  i.e.,  the  capillary  pressure 
driven  flow  is  not  the  only  physical  phenomenon  in  cathode. 

Both  experimental  observations  [4,5]  and  numerical  simulation 
[3]  show  that  the  vapor  condensation  (liquid-vapor  phase  change) 
plays  a  pivotal  role  in  the  water  configuration.  The  vapor  conden¬ 
sation  is  affected  by  the  temperature  difference  in  the  GDL  and 
the  local  wettability.  The  liquid  water  configuration  dominated 
by  phase  change  differs  completely  from  that  dominated  by  the 
capillary  pressure  driven  flow.  If  it  is  completely  dominated  by 
capillary  pressure,  the  liquid  water  will  flow  from  the  higher  sat¬ 
uration  region  (GDL-MPL  interface)  to  the  lower  saturation  region 
(GDL-channel  interface)  because  a  large  capillary  pressure  is  given 
in  the  high  saturation  region.  Continuous  liquid  water  flow  in 
micro-pores,  one  of  the  basic  assumptions  in  capillary  pressure 
driven  flow  in  porous  media,  is  never  observed  in  situ  in  a  PEM  fuel 
cell  system.  On  the  contrary,  a  series  of  discrete  small  droplets  in 
the  GDL  are  often  observed  using  the  environmental  scanning  elec¬ 
tron  microscope  (ESEM)  [6,7].  Although  the  transport  process  in  the 
GDL  cannot  be  clearly  observed  in  situ  in  a  PEM  fuel  cell  system,  it 
can  be  verified  that  the  liquid  water  in  the  GDL  is  discrete  even 
when  the  water  droplets  intrude  into  the  gas  flow  channel  [4,5,8]. 
This  is  inconsistent  with  the  basic  assumption  of  capillary  pressure 
driven  flow.  Therefore,  the  phase  change  may  be  another  key  factor 
to  affect  the  liquid  water  configuration  besides  the  capillary  pres¬ 
sure.  ESEM  imaging  [6,7]  and  synchrotron  X-ray  radiography  [4,5] 
provide  a  high  spatial  resolution  technique  to  investigate  droplet 
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formation  at  the  microscale,  however,  the  ability  to  simulate  fuel 
cell  operating  conditions  has  not  yet  been  reported  in  the  literature 
[9]. 

The  hydrophobic  treatment  technology  [10]  for  the  GDL  fibers 
can  reduce  the  hydrophilic  regions  where  the  vapor  condensation 
easily  occurs.  However,  the  hydrophobic  treatment  increases  the 
hydraulic  pressure  for  driving  the  liquid  water  flow  through  the 
micro-pores.  It  has  been  shown  that  an  appropriate  content  of 
Teflon  (5-30  wt%)  in  the  GDL  can  increase  the  performance  of  PEM 
fuel  cell  [10-12].  There  are  some  other  factors  controlling  the  liquid 
configuration  besides  the  fraction  of  the  hydrophobic  fibers  (e.g., 
fibers  coated  with  Teflon),  for  example,  the  non-uniform  distribu¬ 
tion  of  the  hydrophobic  fibers  discussed  by  Sinha  and  Wang  [13].  It 
is  generally  accepted  that  a  higher  pressure  is  required  to  push  the 
water  through  a  hydrophobic  GDL  than  through  a  hydrophilic  GDL 
[14],  and  the  liquid  water  preferentially  flows  through  the  con¬ 
nected  hydrophilic  pore  network  of  a  mixed-wet  GDL  [13].  This 
paper  will  study  the  quantitative  relationship  between  the  liq¬ 
uid  water  transport  resistance  and  the  pore  properties  in  the  GDL, 
including  the  geometry  parameters  and  wettability. 


2.  Analysis  methods 

The  gas  diffusion  layer  is  generally  made  of  either  a  woven  car¬ 
bon  cloth  or  a  nonwoven  carbon  paper.  It  has  usually  been  modeled 
as  a  homogeneous  porous  medium  with  a  constant  permeability  in 
the  published  literatures  of  PEM  fuel  cells  [15-17].  The  methods 
used  in  those  literatures  show  a  high  numerical  efficiency  since 
they  do  not  require  very  fine  grid  to  capture  the  microstructure 
of  the  GDL.  However  these  macroscopic  models  fail  to  take  into 
account  the  influence  of  the  local  wettability  and  structural  mor¬ 
phology  which  affect  the  cell  performance  [18].  In  fact,  most  of 
the  GDLs  are  not  homogeneous.  For  a  carbon  cloth  GDL,  the  size 
of  the  carbon  fiber  bundle  is  about  200  [Jim,  close  to  the  thick¬ 
ness  of  the  GDL  (~300  pun).  For  a  carbon  paper  GDL,  the  largest 
pore  size  reaches  over  100  pun,  also  close  to  the  thickness  of  the 
GDL  [1,19,20].  According  to  the  experimental  observation  [14],  the 
unsaturated  flow  in  the  GDL  is  affected  by  the  wettability  of  the 
local  fiber  surface,  the  pore  size  distribution,  and  the  fiber  diame¬ 
ter,  etc.  A  simplified  model  is  here  proposed  to  study  the  effects  of 
capillary  pressure,  gravity,  vapor  condensation  speed,  wettability 


and  microstructures  of  the  GDL  on  the  liquid  configuration.  Both 
the  analytical  and  LBM  models  are  adopted  in  the  analysis. 


2.1.  Analytical  method 


The  analytical  analysis  is  based  on  the  Young-Laplace  law, 
which  gives  the  relationship  of  the  pressure  difference  across  the 
liquid-gas  interface,  the  interface  shape  and  surface  tension: 

A  P  =  yK=Z-  (1) 

Kd 

where  AP  is  the  pressure  difference  across  the  fluid  interface,  y  is 
the  surface  tension,  /c  is  the  curvature  of  the  interface,  and  Rd  is  the 
curvature  radius  of  liquid-gas  interface. 

The  difference  of  gas  and  liquid  pressure  within  a  porous 
medium,  i.e.,  the  capillary  pressure,  is  balanced  by  the  tension  of 
the  liquid-gas  interface  in  the  pores.  In  the  previous  studies,  the 
GDL  pores  are  often  idealized  as  the  equivalent  cylindrical  capil¬ 
laries  (usually  called  equivalent  cylindrical  capillary  model)  [21]. 
Based  on  the  Young-Laplace  law,  the  relationship  between  capil¬ 
lary  pressure  and  equivalent  pore  size  is  described  as: 


P  P  _  2ycos(6>) 

Pc-P,-Pg-  — 


(2) 


where  Rp  is  the  pore  radius,  P\  is  the  pressure  of  liquid  phase,  Pg  is 
the  pressure  of  gas  phase,  0  is  the  contact  angle  of  liquid  water  on 
the  fiber  surface. 

Based  on  the  equivalent  cylindrical  capillary  model,  a  zero  cap¬ 
illary  pressure  indicates  that  all  the  pores  are  filled  with  liquid  for 
a  hydrophilic  GDL,  but  all  the  pores  are  drained  for  a  hydropho¬ 
bic  GDL.  However,  to  the  best  of  our  knowledge,  the  phenomenon 
that  liquid  water  is  imbibed  into  the  GDL  without  a  pressure  gradi¬ 
ent  is  never  observed  in  the  experiment.  In  the  present  paper,  we 
proposed  an  equivalent  capillary  model  of  fiber-fence  structure  to 
expound  the  fundamental  mechanisms  of  liquid  water  transport 
in  the  GDL  (see  Fig.  1).  The  continuous  solid  wall  in  the  equivalent 
cylindrical  model  [21  ]  is  here  replaced  by  the  microscale  fiber-fence 
structure.  Since  liquid  water  transports  in  a  quasi-static  process  at 
low  values  of  Capillary  number  and  Reynolds  number,  the  dynamic 
effect  is  not  considered  in  the  analytical  analysis.  According  to  the 
Young-Laplace  law  and  the  geometry  relationship  between  the 


Fig.  1.  Schematic  of  the  equivalent  capillary  model  of  the  fiber-fence  structure. 
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carbon  fibers  and  the  liquid-gas  interface,  the  gas-liquid  interface 
profile  is  obtained  (see  Fig.  1): 


[rf  cos(of)  -  Rp] 
cos(a  +  0) 


(3) 


h  =  Rd[  1  -  sin(a  +  0)]  +  rf  sin(aO  +  nL 


(4) 


where  6  is  the  contact  angle  of  the  liquid  water  on  the  fiber  surface, 
2RP  the  width  of  the  capillary  tube,  rf  the  radius  of  the  fiber,  n  the 
fiber  layer  number  that  the  water  surface  reaches  (equaling  2  in  the 
condition  as  shown  in  Fig.  1 ),  L  the  distance  of  two  neighboring  layer 
fibers,  h  the  center  height  of  the  liquid  level.  The  pore  size  2 Rp  must 
be  greater  than  the  fiber  distance  L  -  2rf,  otherwise,  a  side  leakage 
of  the  liquid  water  may  occur  between  the  two  neighboring  fibers. 

During  the  liquid  surface  imbibition  (rising  as  shown  in  Fig.  1  a),  a 
transition  occurs  when  the  liquid  surface  touches  the  inner-bottom 
surface  of  the  upper  layer  fibers.  But  during  the  liquid  surface 
drainage  process  (going  down  as  shown  in  Fig.  lb),  the  transi¬ 
tion  occurs  when  the  liquid  surface  almost  leaves  the  inner-bottom 
surface  of  the  touching  layer  fibers.  Here  we  defined  a  advancing 
critical  radius  in  the  imbibition  process,  Rd,ac>  equaling  the  cur¬ 
vature  radius  of  the  liquid  water  surface  at  the  transition  point 
(Fig.  la).  During  the  drainage  process,  a  receding  critical  radius, 
ftd.ro  occurs  when  the  contact  line  is  just  separated  from  the  fiber 
surface  (Fig.  lb).  In  these  two  critical  conditions,  the  geometrical 
relationships  are  given  as  follows: 

Pac  -  Rd,ac)  -  (n  +  1  )L}2  +  R2  =  (Rd,ac  +  rf)2  (5) 


(Peng-Robinson)  EOS  [30]  gives  a  better  result  with  low  spurious 
currents. 

In  this  paper,  the  computation  code  of  the  SC  LBM  was  designed 
by  us  and  was  validated  through  careful  comparisons  with  the 
analytical  results  of  the  equivalent  capillary  model  of  fiber-fence 
structure.  Here,  we  implement  the  SC  LBM  in  two-dimension  flow 
analysis  for  a  single-component  multiphase  flow  system.  In  the 
D2Q9  (two-dimension  and  nine  discrete  velocities)  model,  a  stan¬ 
dard  LBM  with  the  BGK  (Bhatnagar-Gross-Krook)  collision  term  is 
adopted,  describing  the  evolution  of  particle  distribution  function 
(PDF )fi  in  site  x  and  time  t  as  follows: 

fi(x  +  etst,  t  +  St)  =  /,(x,  t )  -  L [/f(x,  t)  -/;eq(x,  0], 

i  =  0,1,2,...,  8  (9) 

where  8t  is  the  time  step,  subscript  i  indicates  the  discrete  veloc¬ 
ity  direction,  r  is  the  single  relaxation  time,  having  the  relation 
with  the  kinematic  viscosity:  v  =  (r  -  0.5 )cj8t,  cs  =  c/^/3,  where 
c  =  8xl8t  is  the  ratio  of  lattice  spacing  8x  and  time  step  8t.  is 
equilibrium  PDF,  and  can  be  calculated  as 

'  erueq  (e,-  ueq)2  (ueq)2" 
cj  2c f  2 cf  _  ’ 

r  4/9,  i  =  0 

&>,-  =  <j  1/9,  i  =  l-4  (10) 

^  1/36,  i  =  5-8 


t)  =  a>ip 


drc  +  Prc  —  ~2  (6) 

where  the  subscripts  ac  and  rc  denote  the  advancing  critical  value 
and  receding  critical  value  respectively. 

From  Eq.  (3),  one  obtains: 

dRd  _  rf  sin(0)  -  Rp  sin(a  +  0) 
da  cos2{a  +  0) 

When  dRd/da  =  0,  i.e., 


In  Eqs.  (9)  and  (10),  e,-  is  the  discrete  velocities,  given  by 


0  10-1  0  1-1-1  1 
0  0  10-111-1-1 


(11) 


The  macroscopic  density  and  the  modified  macroscopic  velocity 
are  given  as  follows: 

P  =  YJi  (12) 


rf  sin(0)  -  Rp  sin{a  +  0)  =  0  (8) 

combining  Eqs.  (2)  and  (3),  the  minimum  curvature  radius  and  the 
corresponding  maximum  capillary  pressure  are  obtained.  The  max¬ 
imum  capillary  pressure  is  just  the  hydraulic  pressure  required  for 
the  liquid  water  to  penetrate  into  the  GDL. 


ueq  =YJrei  +  r1St1F_ 
P 


(13) 


In  the  single-component  multiphase  LB  model  proposed  by  Shan 
and  Chen,  a  simple  long-range  interaction  force  between  the  par¬ 
ticle  at  site  x  and  the  particle  at  site  x'  is  introduced  as: 


2.2.  Multiphase  lattice  Boltzmann  method 

In  the  above  section,  we  propose  an  analytical  method  to  ana¬ 
lyze  the  transport  of  liquid  water  in  the  GDL  induced  only  by  the 
capillary  pressure.  In  the  present  section,  we  will  use  the  lattice 
Boltzmann  method  (LBM)  [22],  to  study  the  transport  of  liquid 
water  in  the  GDL  when  the  vapor  condensation  and  stochastic 
microstructure  are  considered.  The  LBM  has  been  accepted  as  a 
new  computational  tool  for  two-phase  flow  with  phase  change  in 
complex  geometries,  and  has  been  used  for  simulating  the  liquid 
water  in  the  GDL  in  order  to  obtain  a  fundamental  understanding 
[23-27].  Here  we  use  the  multiphase  and  single-component  LBM 
originally  proposed  by  Shan  and  Chen  (SC)  [28]  to  simulate  the 
water  flow  in  the  porous  GDL.  The  SC  multiphase  LBM  used  a  con¬ 
cept  of  the  microscopic  interactions  between  the  fluid  particles  at 
the  neighboring  lattices  by  adding  the  additional  forcing  term  to 
the  velocity  field.  The  interaction  force  is  controlled  by  the  equa¬ 
tion  of  state  (EOS).  The  spurious  currents  (the  unphysical  velocities 
generated  in  the  simulation)  are  ubiquitous  in  the  SC  model  and 
cannot  be  separated  from  the  real  flow  velocities.  Yuan  and  Schae¬ 
fer  [29]  discussed  some  choices  of  EOS  and  proposed  that  the  P-R 


F(x)  = 


Vf(x)Ex'®iVf(x')(x'  -  X) 


(14) 


where  is  called  the  “effective  mass”  or  “pseudo  potential  func¬ 
tion”  and  is  defined  as  a  function  of  x  as: 


x//(x)  =  a/2 (p(x)cs2  -P(/0(x))) 


(15) 


The  pressure  p  of  non-ideal  gas  is  described  by  the  EOS. 

The  single-component  multiphase  LB  model  is  often  limited  to 
a  small  density  ratio  because  numerical  instability  may  appear  in 
case  of  a  large  density  ratio.  Yuan  and  Schaefer  [29]  discussed  some 
choices  of  EOS  and  proposed  that  the  P-R  EOS  gives  better  results 
in  cases  of  large  density  ratio.  The  pressure  of  non-ideal  gas  is 
described  by  the  P-R  EOS  as: 


_  pRT  aa{T)p2 

P=  \-bp~  1  +2b/o-b2/02 


(16) 


«(T)  = 


1  +(0.37464  +  1.54226« -0.26992<w2)x  (l  -  ^ T/T( 


-i  2 


(17) 
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where  a  =  0.445724 R2Tc2/pc,  b  =  0.08778 RTc/pc,  Tc  is  the  critical 
temperature,  pc  is  the  critical  pressure,  R=  1  in  the  present  simula¬ 
tion,  co  is  the  acentric  factor  equaling  0.344  for  water  [29]. 

The  wetting  phenomena  of  the  liquid-solid  interface  can  be 
implemented  by  appointing  a  pseudo  effective  mass  ifsw  to  simulate 
the  interaction  force  (similar  to  the  fluid-fluid  interaction  [31])  at 
the  fluid-solid  boundary.  The  wetting  property  of  the  solid  bound¬ 
ary  is  given  by  adjusting  the  ifsw  from  i//g  (effective  mass  of  gas 
phase)  to  i/q  (effective  mass  of  liquid  phase). 

The  analysis  region  is  shown  in  Fig.  2.  Similar  to  the  homogenous 
macroscopic  2D  model  [32,33],  it  contains  the  flow  channel  and 
GDL.  Boundary  B4  is  treated  as  a  symmetric  boundary,  B1  is  the 
interface  of  the  GDL  and  MPL,  B2  is  the  surface  of  the  bipolar  plate, 
and  B3  is  the  surface  of  the  fibers.  The  MPL  is  treated  as  a  boundary 
in  the  model.  No-slip  boundary  condition  is  used  at  Bl,  B2,  and  B3. 
The  diameter  of  the  carbon  fiber  is  1 0  p,m,  and  the  porosity  is  ~80%, 
which  is  close  to  the  real  porosity  of  Toray  TGP-H  carbon  paper. 
The  width  of  the  analysis  region  is  1  mm,  and  thickness  of  the  GDL 
is  150  |jim.  The  lattice  spacing  Sx  equals  0.5  pum.  At  the  initial  step, 
the  carbon  fibers  are  placed  either  randomly  or  regularly,  while  the 
whole  region  is  occupied  by  the  saturated  vapor.  Due  to  the  water- 
free  spots  caused  by  the  Teflon,  the  water  in  the  MPL  appears  in  the 
form  of  vapor.  Source  terms  (denoting  the  vapor  diffused  into  the 
GDL  from  MPL)  are  added  to  the  evolution  equations  of  the  lattices 
at  Bl.  The  water  content  increases  continuously  with  a  very  small 
velocity.  When  the  liquid  water  droplet  intrudes  into  the  channel 
and  becomes  big  enough,  it  is  blown  away  along  the  channel.  In  the 
simulation  process  in  the  simplified  2D  model,  the  water  droplet 
in  the  channel  is  cleared  away  instantaneously.  However,  how  the 
droplet  to  move  away  after  its  formation  in  the  channel  is  out  of 
the  scope  of  the  present  paper. 

In  the  LBM  analysis,  the  water  contact  angle  of  the  bipolar  plate 
is  assumed  to  be  80°  [34],  and  the  contact  angle  of  the  MPL  surface 
(Bl  in  Fig.  2)  is  assumed  to  be  150°  [35].  As  for  the  wettability  of  the 
fibers  in  the  GDL,  a  contact  angle  of  1 10°  corresponds  to  the  carbon 
fibers  coated  by  Teflon,  but  a  contact  angle  of  80°  corresponds  to 
the  untreated  carbon  fibers  [36]. 

We  use  a  constant  pore  size  distribution  model  (i.e.,  the  GDL 
consists  of  uniformly  arrayed  fibers)  to  analyze  the  effect  of  the 
mixed-wettability,  but  use  the  random  pore  size  distribution  model 
to  analyze  the  process  of  liquid  water  removal  from  the  MPL  to  the 
flow  channel.  The  effect  of  pore  size  distribution  is  here  discussed 
qualitatively,  such  as  regular  distribution  and  random  distribution. 

3.  Results  and  discussion 

3.1.  Capillary  pressure  driven  flow  in  porous  media 

Based  on  the  equivalent  capillary  model  of  fiber-fence  structure, 
the  curvature  radius  of  the  liquid-gas  interface  is  determined  from 
Eqs.  (3)  and  (4).  For  each  calculation  point,  the  transition  condition 
is  judged  by  Eq.  (5).  When  Eq.  (5)  is  satisfied,  the  variable  n  will  be 
increased  by  1  in  Eq.  (4).  The  volume  of  the  liquid  water  is  assumed 
unchanged  before  and  after  jumping  of  the  water  surface.  The  cur- 


Fig.  3.  The  liquid  surface  property  vs.  the  height  of  the  liquid  surface  for  different 
wettabilities  of  the  fiber  surfaces  in  the  imbibition  process,  (a)  The  curvature  radius 
and  (b)  the  capillary  pressure.  1  =  20  p,m  and  Rp  =  15  p,m.  The  dashed  line  expresses 
the  transition  process  in  which  the  liquid-gas  interface  touches  the  inner-bottom 
surface  of  upper  layer  of  the  fibers,  n  denotes  the  fiber  layer  number  that  the  water 
surface  reaches. 


vature  radii  of  the  liquid-gas  interface  are  shown  in  Fig.  3a  when 
the  liquid  water  gradually  intrudes  into  the  upper  fiber  layers  in  the 
imbibition  process.  The  change  of  the  curvature  radius  of  the  liquid 
surface  shows  approximately  a  periodical  behavior  except  at  the 
regions  nearby  the  transition  points.  When  n  is  less  than  3  and  the 
fibers  are  hydrophilic,  the  curvature  radius  will  become  a  negative 
value,  i.e.,  the  liquid-gas  interface  is  concave.  However,  the  curva¬ 
ture  is  always  positive  when  the  fibers  are  hydrophobic.  Except  at 
the  regions  nearby  the  transition  points,  the  influence  of  the  fiber 
wettability  on  the  curvature  radius  is  not  as  large  as  expected. 

The  corresponding  capillary  pressures  predicted  by  Eq.  (1)  are 
shown  in  Fig.  3b.  The  curves  show  a  periodic  behavior  and  there 
is  a  maximum  capillary  pressure  in  each  period.  To  drive  the  liq¬ 
uid  water  going  through  the  capillary  pore,  a  hydraulic  pressure 
should  be  applied  to  overcome  the  maximum  capillary  pressure. 
When  the  contact  angle  of  the  fiber  surface  is  changed  from  80°  to 
110°,  the  maximum  capillary  pressure  is  increased  from  4.79  kPa 
to  5.74 kPa.  However,  based  on  Eq.  (2),  the  equivalent  cylindri¬ 
cal  model  gives  the  corresponding  equivalent  contact  angles  of 
120°  and  127°  respectively.  Therefore,  all  of  these  pores  belong  to 
hydrophobic  pores  according  to  the  equivalent  cylindrical  model. 
The  maximum  capillary  pressure  is  always  positive  even  when  the 
contact  angle  is  60°  as  shown  Fig.  4.  The  maximum  capillary  pres¬ 
sure  is  greater  than  1 .7  kPa  when  the  width  (2Rd  as  shown  in  Fig.  1 ) 
of  the  capillary  pore  ranges  from  20  p,m  to  80  p,m  and  the  contact 
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Fig.  4.  The  maximum  capillary  pressure  to  drive  the  liquid  water  flowing  through 
the  pores.  1  =  20  p,m. 


angle  ranges  from  60°  to  1 1 0°.  Therefore,  the  liquid  water  would  not 
intrude  into  the  GDL  without  an  applied  hydraulic  pressure  if  the 
phase  change  and  the  gravity  are  both  negligible.  In  fact,  the  max¬ 
imum  capillary  pressure  is  always  positive.  Now  let  us  discuss  the 
limiting  situation.  For  the  completely  wetting  (zero  contact  angle) 
fibers,  the  water  surface  curvature  is  zero  when  the  height  of  the 
liquid  surface  equals  nl  +  rf.  Then,  the  curvature  is  positive  with  the 
increase  in  the  height  of  the  liquid  surface  before  the  liquid  surface 
touches  the  next  layer  of  fibers  (i.e.,  a  surface  jumping  occurs).  The 
jumping  phenomenon  in  the  capillary  pressure  is  speculated  as  a 
special  phenomenon  occurring  in  the  fiber  porous  media. 

The  equivalent  capillary  model  of  fiber-fence  structure  can  be 
used  to  reveal  the  mechanism  of  the  liquid  water  transport  in  the 
GDL  induced  by  capillary  pressure.  But  effects  of  the  phase  change 
and  the  body  force  cannot  be  included  in  the  analytical  solution. 
Therefore,  a  numerical  solution  for  the  LBM  will  be  carried  out 
in  the  following  analysis  to  study  the  effects  of  gravity  and  phase 
change,  and  the  transport  process  of  water.  The  analytical  and  LBM 
solutions  based  on  the  equivalent  capillary  model  of  fiber-fence 
structure  are  compared  in  Fig.  5a  and  b  for  the  different  fiber  wet¬ 
tabilities.  In  Fig.  5a  we  give  the  comparison  between  the  analytical 
solution  and  LBM  solution  for  a  large  range  of  water  level  height 
(5  layers  of  fibers)  when  the  contact  angle  0  =  80°.  In  Fig.  5b  we 
give  a  detailed  comparison  of  the  analytical  solutions  and  the  LBM 
solutions  for  0  =  60°,  80°,  110°  and  140°  respectively  at  the  first  two 
layers  of  fibers.  Good  agreements  occur  between  the  LBM  simula¬ 
tions  and  the  analytical  solutions  except  at  the  regions  nearby  the 
transition  points.  The  LBM  simulation  shows  a  quasi-static  process 
of  the  liquid  surface  elevation,  but  the  analytical  solution  gives  a 
completely  static  process. 


3.2.  Gravitational  effect 


The  Bond  number  (denoted  here  by  B0 ),  a  dimensionless  number 
expressing  the  ratio  of  the  body  forces  to  the  surface  tension  force, 
is  expressed  as: 


(18) 


where  p  is  the  density  of  the  liquid,  g  the  gravity  acceleration,  l  the 
characteristic  length(the  radius  of  a  droplet  or  the  radius  of  a  capil¬ 
lary  tube),  y  the  surface  tension  of  the  interface.  When  the  droplet 
radius  is  100|jim  in  the  GDL,  the  Bond  number  is  about  0.0015. 
Therefore,  effect  of  the  gravity  on  the  capillary  flow  in  the  GDL  is 
very  small  and  is  not  considered  in  the  present  paper.  Flowever, 


Fig.  5.  Comparison  between  the  analytical  solution  and  the  LBM  simulation  for  the 
surface  curvature  radius  of  the  liquid  water  for  different  contact  angles  as  indicated 
in  the  figures  (L  =  20  |jim  and  Rp  =  15  |xm):  (a)  the  analytical  solution  and  LBM  method 
for  contact  angle  0  =  80°  for  a  large  range  of  water  level  height  and  (b)  the  detailed 
comparison  of  the  analytical  solution  and  the  LBM  solutions  for  6  =  60°,  80°,  110° 
and  140°  respectively  at  the  first  two  layers  of  fibers. 


the  largest  droplet  intruding  in  the  gas  flow  channel  will  be  greater 
than  1  mm,  and  the  corresponding  Bond  number  is  greater  than 
0.136.  In  this  case,  the  gravity  may  affect  the  water  droplet  config¬ 
uration  and  detachment.  The  performance  of  the  PEM  fuel  cell  is 
gravity-dependent  as  observed  by  Kimball  et  al.  [37]. 

3.3.  Effect  of  vapor  condensation 

The  water  removal  induced  by  vapor  condensation  (or  phase 
change)  is  affected  by  the  wettability  of  GDL  fibers,  the  vapor  diffu¬ 
sion  velocity,  and  the  temperature  gradient.  The  temperature  effect 
is  not  considered  in  the  present  study  due  to  the  limitation  of  SC 
LBM  despite  a  little  temperature  difference  (~3  K)  may  have  obvi¬ 
ous  effect  on  the  water  distribution  as  reported  in  Ref.  [3].  Up  to 
now,  the  single-component  multiphase  SC  LBM  still  cannot  simu¬ 
late  the  situation  at  low  temperature  (<373 1<)  and  high  density  ratio 
due  to  the  numerical  instability  at  large  density  ratio.  The  temper¬ 
ature  studied  in  the  present  paper  is  0.8TC  (TC  =  647I<).  Therefore, 
the  analysis  for  the  phase  change  process  is  qualitative. 

When  the  vapor  condensation  effect  is  considered  in  the  LBM 
model,  the  liquid  water  distribution  is  dominated  by  the  vapor  dif¬ 
fusion  velocity  and  the  water  production  rate.  Fig.  6  shows  the 
liquid  water  distribution  for  several  values  of  the  source  terms 
when  the  average  liquid  water  saturation  in  the  GDL  is  ~15%.  The 
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Fig.  6.  The  liquid  water  configuration  for  different  source  terms:  (a)  0.005,  (b)  0.002, 
(c)  0.001,  and  (d)  0.0005.  All  the  fibers  are  hydrophilic  and  give  a  random  distribu¬ 
tion.  The  blue  regions  denote  the  liquid  water  droplets.  (For  interpretation  of  the 
references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of 
the  article.) 

source  terms  are  0.005,  0.002,  0.001,  and  0.0005  respectively  for 
each  time  step.  A  large  velocity  of  phase  change  occurs  where  the 
vapor  pressure  is  high.  Due  to  the  flow  and  diffusion  resistances, 
the  vapor  pressure  near  the  MPL-GDL  interface  (B1 )  is  higher  than 
that  at  the  top  surface  of  the  GDL,  and  therefore  the  vapor  will 
condense  at  the  small  pores  near  Bl.  Generally,  the  small  pores 
are  surrounded  by  more  carbon  fibers  than  the  large  pores,  and 
thus  give  more  hydrophilic  property  than  the  large  pores.  When 
the  source  term  is  reduced,  the  pressure  difference  between  the 
top  and  bottom  regions  will  decrease  due  to  the  abundant  time  for 
vapor  diffusion.  Therefore,  the  condensation  occurs  at  the  whole 
GDL  region  but  not  only  at  the  region  near  Bl.  Based  on  the  LBM 
simulations  for  different  source  terms,  we  found  that  the  two  main 
positions  of  liquid  water  condensation  are  located  close  to  the 
interface  of  the  rib  land  and  the  MPL-GDL  interface  (Bl).  This  is 
supported  by  the  experimental  observations  [4,5]  and  numerical 
simulation  [3].  Even  though  the  temperature  effect  is  not  consid¬ 
ered,  the  present  simulation  shows  that  a  hydrophilic  land  surface 
can  easily  induce  the  vapor  condensation  occurring  under  the  land. 
This  agrees  with  the  experimental  observations  [4,5]  although  the 
physical  reasons  in  the  experiments  may  contain  the  temperature 
effect.  Due  to  little  physical  basis  of  the  phase  change  velocity  in  the 
SC  LBM  model,  the  source  term  is  still  difficult  to  be  translated  in  to 
real  world  units.  The  current  density  is  3.72  x  10-5/<$t  (Ccm-2) 
when  the  value  of  the  source  term  is  0.001  (If  the  time  step  is 
0.5  pus,  the  current  density  is  74.4  A  cm-2.).  The  real  time  step  is 


still  difficult  to  be  given  due  to  unreal  phase  change  velocity  in  the 
model. 

3.4.  Liquid  water  formation  and  transport 

What  is  discussed  above  shows  that  the  capillary  driven  flow 
and  phase  change  are  the  two  main  physical  principles  of  the  water 
transport  in  GDL.  Here  we  will  analyze  the  transport  process  of  the 
liquid  water  from  MPL  towards  the  flow  channel  and  the  effect 
of  the  spatial  mixed-wettability  distribution  on  the  liquid  water 
configuration. 

Fig.  7(a)-(f)  shows  the  liquid  water  configuration  when  the 
carbon  fibers  distribute  regularly  but  have  different  surface  wetta¬ 
bilities.  The  average  saturation  in  the  GDL  is  45%,  and  the  porosity  is 
~80%.  The  insets  at  the  top  right  corners  (lands)  in  the  figures  show 
the  fiber  distribution  types:  (a)  all  the  fibers  are  hydrophilic;  (b)  all 
the  fibers  are  hydrophobic;  (c)-(f)  the  hydrophobic  fiber  fraction 
is  taken  as  a  constant  (50%),  but  the  hydrophobic  and  hydrophilic 
regions  are  different.  The  black  dots  denote  the  hydrophilic  fibers, 
while  the  orange  dots  denote  the  hydrophobic  fibers.  In  such  a  case, 
we  found  that  the  liquid  water  distribution  is  completely  different. 
Therefore,  the  fraction  of  hydrophobic  fibers  cannot  itself  describe 
the  transport  ability  of  the  liquid  water  through  the  GDL.  The  dif¬ 
ference  of  the  maximum  pressure  between  the  hydrophilic  and 
hydrophobic  pores  induces  the  liquid  water  preferentially  flowing 
into  the  hydrophilic  pores.  The  hydrophobic  fibers  surrounding  the 
hydrophilic  regions  give  a  large  resistance  to  the  liquid  water  flow. 
The  fiber  distribution  patterns  shown  in  Fig.  7(d)  and  (f)  seem  more 
efficient  for  the  water  management  and  the  cell  efficiency  since  the 
gas  paths  in  the  GDL  are  not  hindered  by  the  liquid  water.  The  liq¬ 
uid  water  easily  congregates  at  the  land-GDL  interface  due  to  the 
hydrophilic  surface  of  the  land  (see  Fig.  7(a)-(d)).  The  aforemen¬ 
tioned  analysis  is  carried  out  based  on  the  assumption  of  regularly 
distributed  carbon  fibers. 

Fig.  8(a)-(d)  shows  the  liquid  water  configuration  for  different 
average  saturations.  In  Fig.  8  all  the  fibers  are  hydrophilic.  The  pro¬ 
cess  of  liquid  water  transport  from  the  MPL  to  the  GDL  can  thus  be 
summarized  as  follows:  initially,  the  liquid  water  condenses  at  the 
small  pores,  then  some  of  the  water  droplets  grow  up,  but  others 
may  either  stop  growing  up  or  go  off.  Once  they  grow  up  to  big 
enough  to  touch  more  hydrophilic  fibers,  they  will  never  vanish. 
After  they  are  surrounded  by  too  many  fibers  the  droplets  will  stop 


Fig.  7.  The  liquid  water  configuration  for  different  fiber  distributions  and  wettabilities.  The  average  saturation  is  45%  and  the  porosity  is  ~80%.  The  square  insets  at  the  top 
right  corners  show  the  distribution  types  of  the  hydrophobic  carbon  fibers  (coated  with  Teflon).  The  black  dots  denote  the  hydrophilic  fibers,  while  the  orange  dots  denote 
the  hydrophobic  fibers.  The  blue  regions  denote  the  liquid  water.  The  land  surfaces  from  (a)  to  (d)  are  hydrophilic,  while  the  land  surfaces  in  (e)  and  (f)  are  hydrophobic.  (For 
interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 
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Fig.  8.  The  liquid  water  configuration  and  the  evolution  process  for  different  average  saturations:  (a)  8%,  (b)  15%,  (c)  39%,  and  (d)  47%.  The  blue  regions  denote  the  liquid 
water  droplets.  All  the  fibers  are  hydrophilic.  The  red  dots  in  (c)  denote  the  fibers  starting  to  form  the  initial  condensation.  (For  interpretation  of  the  references  to  color  in 
this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 


growing  up,  such  as  the  droplet  at  the  bottom  right  corner.  This  is 
because  that  a  high  hydraulic  pressure  is  required  for  liquid  water 
to  flow  through  the  small  pores.  Another  interesting  phenomenon 
is  that  when  the  liquid  water  flows  through  the  GDL  into  the  flow 
channel,  the  nearby  droplets  will  stop  growing  up.  As  shown  in 
Fig.  8(c)  and  (d),  all  the  droplets  under  the  flow  channel  stop  grow¬ 
ing  up  when  the  liquid  water  intrudes  into  the  flow  channel.  Most 
of  the  water  droplets  will  encounter  the  dead  ends  in  the  GDL  and 
occupy  more  space,  especially  at  the  regions  under  the  land.  This 
agrees  with  the  pore-network  modeling  simulation  of  Sinha  and 
Wang  [38]. 

Fig.  8(d)  shows  the  process  of  the  movement  paths  of  the  liq¬ 
uid  water.  21  red  points  (fibers  with  red  color)  indicate  the  initial 
positions  of  condensation  occurring.  But  only  10  of  them  grow  up. 
The  lines  and  arrowheads  show  the  routes  and  directions  of  the 
liquid  water  flow.  The  liquid  water  will  flow  through  the  big  pores 
between  fibers.  This  is  because  the  curvature  of  the  droplet  at  the 
big  pore  is  smaller  than  that  at  the  small  pore,  and  therefore  the 
pressure  difference  for  liquid-gas  interface  movement  in  the  big 
pore  is  also  smaller  than  that  in  the  small  pore.  Therefore,  for  a 
good  design  of  GDL,  the  following  two  conditions  are  required:  (a) 
there  are  some  connected  big  pores  from  MPL  to  GDL,  and  (b)  there 
are  enough  small  hydrophobic  pores  to  afford  space  for  the  reactant 
gas  to  flow. 

The  temperature  gradient  is  neglected  in  the  present  analysis. 
This  simplification  induces  an  underestimation  of  the  velocities 
of  vapor  condensation  in  the  channel  and  the  region  nearby  the 
land.  If  the  temperature  gradient  exists  in  the  cathode,  the  droplet 
would  congregate  where  the  temperature  is  lower  for  the  same 
wettabilities  of  fiber  surfaces. 

4.  Conclusions 

Understanding  the  process  of  the  liquid  water  droplet  formation 
and  transport  in  the  GDL  is  a  key  to  unraveling  the  origin  and  devel¬ 
opment  of  flooding.  In  the  present  paper,  the  equivalent  capillary 
model  of  fiber-fence  structure  and  lattice  Boltzmann  method  are 
applied  to  simulate  the  micro-scale  multiphase  flow  in  the  GDL. 
The  possible  contributors  to  liquid  water  configuration  are  dis¬ 
cussed,  such  as  capillary  pressure,  gravity  and  phase  change.  We 
also  qualitatively  discuss  the  effect  of  fiber  distribution  and  spatial 
mixed-wettability  distribution  on  the  liquid  water  transport  in  the 
microstructure  of  the  GDL.  We  delineate  the  fundamental  process 
of  the  liquid  water  transport  in  the  GDL.  The  main  findings  and 
conclusions  are  as  follows: 

(1)  A  positive  hydraulic  pressure  is  required  for  driving  the  liq¬ 
uid  water  to  flow  through  the  pores  even  for  the  large  pores 
surrounded  by  hydrophilic  fibers. 


(2)  The  pores  under  the  rib  land  are  easily  blocked  by  liquid  water 
due  to  the  hydrophilic  surface  of  the  land  and  the  high  flow 
resistance  of  the  liquid  water  moving  into  the  flow  channel. 

(3)  The  fraction  of  the  hydrophobic  fibers  itself  cannot  describe  the 
transport  ability  of  the  liquid  water  through  the  GDL.  It  is  found 
that  the  hydrophobic  treatment  pattern  plays  an  important  role 
in  the  liquid  water  configuration. 

(4)  After  the  liquid  water  flows  through  the  GDL  into  the  flow  chan¬ 
nel,  the  nearby  droplets  will  almost  stop  grow  up,  where  a 
gaseous  flow  path  will  be  supplied. 
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